home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / equal_variance.pro < prev    next >
Text File  |  1997-07-08  |  10KB  |  352 lines
  1. ; $Id: equal_variance.pro,v 1.3 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7.  pro equal_variance, X1, FC, DF, P, One_Way = One,  $
  8.                      Unequal_One_Way = Unequal,$
  9.                      Two_Way = Two, Interactions_Two_Way  = $
  10.                      Interactions, Hartley = Hart, Bartlett = Bart, $
  11.                      Box = Bx, Levene = Lev, Missing =M,      $
  12.                      Group_No = gn, Block = B,  TInteraction=IN
  13. ;+
  14. ; NAME:
  15. ;    EQUAL_VARIANCE
  16. ;
  17. ; PURPOSE: 
  18. ;    To test treatments, blocks and interactions, as applicable, for 
  19. ;    variance equality before performing an anova test.
  20. ; 
  21. ; CATEGORY:
  22. ;    Statisitics.
  23. ;
  24. ; CALLING SEQUENCE:
  25. ;    EQUAL_VARIANCE, X [, FC, DF, P]
  26. ;
  27. ; INPUTS:  
  28. ;    X:    two or three dimensional array of experimental data.
  29. ;
  30. ; KEYWORDS: 
  31. ;    BLOCK:    a flag, if set by user, to test homogeneity of block 
  32. ;        variances. Default is treatment variances.  Block should 
  33. ;        only be set when either the keyword TWO_WAY or 
  34. ;        INTERACTIONS_TWO_WAY is set.
  35. ;
  36. ; TINTERACTION:    a flag, if set by the user, to test homogeneity of variance
  37. ;        for all treatment/block combinations.
  38. ;
  39. ; Type of design structure. Options are
  40. ;      ONE_WAY:    if set,  1-way anova.
  41. ;
  42. ; UNEQUAL_ONE_WAY:  if set,  1-way ANOVA with unequal sample sizes.
  43. ;
  44. ; TWO_WAY:    if set,  2-way ANOVA.
  45. ;
  46. ; INTERACTIONS_TWO_WAY:  if set,  2-way ANOVA with interactions.
  47. ;
  48. ; One and only one of the above options may be set.
  49. ; Type of test to be used to test variance equality.
  50. ; Options are:
  51. ;              
  52. ;      HARTLEY:    Hartley's F-Max test. Sample sizes should be equal and data
  53. ;        nearly normally distributed.
  54. ;
  55. ;     BARTLETT:    Bartlett's test for nearly normally distributed data and not
  56. ;        necessarily equal sample sizes.
  57. ;
  58. ;    BOX:    Box's test. Robust for large data sets, data not necessarily
  59. ;        normally distributed.
  60. ;    
  61. ;    LEVENE:    Levene's test.              
  62. ;
  63. ; One and only one test may be selected.
  64. ;
  65. ;      MISSING:    place holder for missing data or undefined if no data is 
  66. ;        missing.
  67. ;
  68. ;     GROUP_NO:    number of groups to use in Box test.
  69. ;
  70. ; OUTPUT:     
  71. ;    FC:    value of statistic computed by each test.  If Bartlett's is 
  72. ;        selected, FC has chi square distribution.  Otherwise, FC has 
  73. ;        F distribution.
  74. ;
  75. ;    DF:    degrees of freedom. 
  76. ;        Hartley: DF = [max sample size -1, number of treatments],
  77. ;        Bartlett: DF = number of treatments,
  78. ;        Box,Levene: DF =[numertor DF,denominator DF]
  79. ;
  80. ;    P:    probability of FC or something greater for Barlett's, Box's 
  81. ;        and Levene's Test.  P = -1 for Hartley's.
  82. ;
  83. ; PROCEDURE:
  84. ;    Milliken and Johnson, "Analysis of Messy Data", Van Nostrand, 
  85. ;    Reinhold, 1984, pp. 18-22.
  86. ;-
  87.  On_Error,2
  88.  
  89.  SX = size(X1)
  90.  
  91.           ;Ascertain design structure.
  92.  ONE_WAY = 0 & ONE_WAY_UNEQUAL=1 & TWO_WAY = 2 & TWO_WAY_INTERACTIONS =3
  93.  keyset = [keyword_set(One), keyword_set(Unequal), $
  94.            keyword_set(Two),keyword_set(Interactions)]
  95.  T  = where(keyset,nkey)     
  96.  
  97.  if nkey NE 1 THEN BEGIN
  98.    print,'equal_variance - must specify one and only one type'
  99.    print, '                of design structure'
  100.   return
  101.  ENDIF
  102.  
  103.  T = T(0)               ;T = selected design structure
  104.  
  105.  if((T EQ TWO_WAY_INTERACTIONS) AND  (SX(0) NE 3) OR $
  106.      ( T NE TWO_WAY_INTERACTIONS AND (SX(0) NE 2))) THEN BEGIN
  107.       print, 'Anova- wrong number of dimensions' 
  108.       return
  109.  ENDIF
  110.  
  111.  
  112.          ; Ascertain test type
  113.  Hartley = 0 & Bartlett = 1 & Box = 2 & Levene = 3
  114.  keyset = [keyword_set(Hart), keyword_set(Bart), $
  115.            keyword_set(Bx),keyword_set(Lev)]
  116.  VT  = where(keyset,nkey)
  117.  
  118.  if nkey NE 1 THEN $
  119.     print, $
  120.       'equal_variance - must specify one and only test type'
  121.  
  122.  VT = VT(0)             ; VT = selected test type
  123.  
  124.  X = float(X1)
  125.  
  126.  SX=Size(X)        
  127.  C=SX(1)      ; Compute number of columns       
  128.  D= 1         ; Default # of replications
  129.  
  130.  
  131.  
  132.  
  133.  if( T EQ TWO_WAY_INTERACTIONS) THEN R = SX(3) else R = SX(2)
  134.  
  135.                                  ;Compute number of rows
  136.  
  137.  
  138.  if R LT 2 OR C LT 2 THEN BEGIN
  139.    print,'equal_variance- data array needs more rows or columns'
  140.    return
  141.  ENDIF
  142.  
  143.  if KEYWORD_SET(B) THEN $
  144.   if T EQ TWO_WAY THEN BEGIN
  145.     X = transpose(X)
  146.     C = R
  147.     R = SX(1)
  148.  ENDIF ELSE if T EQ ONE_WAY or T EQ ONE_WAY_UNEQUAL THEN BEGIN
  149.     print, 'equal_variance- Block test only done with'
  150.     print, '                two_way design structures'
  151.     return
  152.  ENDIF
  153.  
  154. ; Pre-process three dimensional array by converting to a
  155. ; two dimensional array with R*C columns and D rows.
  156.  
  157.  
  158.  if(T EQ TWO_WAY_INTERACTIONS) THEN BEGIN    
  159.  
  160.      D = SX(2)
  161.      if KEYWORD_SET(IN) THEN BEGIN
  162.       X = Fltarr(R*C,D)
  163.       for i = 0L,R*C-1  DO X( i,*) = X1(i mod C,*,i / C) 
  164.       X = float(x)
  165.       C=R*C
  166.  
  167.       R=D
  168.     ENDIF ELSE BEGIN 
  169.         T = ONE_WAY
  170.         if KEYWORD_SET(B) THEN BEGIN
  171.           X = reform(X,C*D,R)
  172.           X = transpose(X)
  173.           temp = C*D
  174.           C = R
  175.           R = temp
  176.         ENDIF ELSE BEGIN
  177.           X = reform(X,C,R*D)
  178.           R = R*D
  179.           D = 1
  180.         ENDELSE
  181.    ENDELSE
  182.  ENDIF
  183.  
  184.  
  185. if (N_Elements(M) NE 0) THEN BEGIN
  186.    if (T EQ TWO_WAY) THEN BEGIN
  187.        print,            $
  188.         "equal_variance- two_way anova without interactions"
  189.        print,"                can't have missing data."
  190.        return
  191.     ENDIF
  192.  
  193.     Pairwise,X,M,YR,YC,NotGood, good               
  194.                                  ; replace missing data by 0,return
  195.                                  ;vectors YR,YC of row,column sizes
  196.      if(N_Elements(NotGood) EQ 0) THEN M=0
  197.        CN = where(YC NE 0,NYC)
  198.        RN = where(YR NE 0,NYR)
  199.        CN1 = where(YC lt 2, NYC1)
  200.        RN1 = where(YR lt 2,NYR1)
  201.  
  202.      if NYC LT 2  OR $
  203.         NYC1  NE 0 THEN BEGIN
  204.        print,"equal_variances, too many missing entries"
  205.       return 
  206.     ENDIF
  207.   ENDIF ELSE M=0
  208.    
  209.  if (VT EQ BARTLETT OR VT EQ HARTLEY ) THEN BEGIN    
  210.                       
  211.     ;compute column variances
  212.  if (M EQ 0) THEN BEGIN
  213.     Mean = (X # replicate(1.0, R))/R
  214.     Var = (X -  MEAN # replicate(1.0, R))^2 
  215.     Var =  Var # replicate(1.0,R)
  216.     Var = Var / (R-1)
  217.  
  218. ;    if (T EQ TWO_WAY) THEN BEGIN       ;add on row variances 
  219. ;       Mean = (replicate(1.0, C) # X)/C
  220. ;       RVAR = (X -  replicate(1.0, C) # MEAN)^2 
  221. ;       RVAR = Replicate(1.0,C) # RVAR
  222. ;      RVar = RVar /(C-1)
  223. ;       Var = [Var, transpose(RVar)]
  224. ;      ENDIF
  225.  
  226.  ENDIF ELSE BEGIN
  227.  
  228.         Mean =  (X # Replicate(1.,R)) /YC  
  229.         Mean = MEAN # replicate(1.0, R)
  230.         VAR = fltarr(C,R)
  231.         VAR(good) = (X(good) - Mean(good))^2
  232.         VAR = VAR # Replicate(1.0,R)
  233.         YC  = YC - 1
  234.         Var = Var/YC
  235.  
  236.      ENDELSE
  237.  
  238.   here = where(VAR eq 0, count)
  239.   if count ne 0 THEN BEGIN
  240.      print, 'equal_variance- cannot handle populations with '
  241.      print, '                0 variance'
  242.      return
  243.   ENDIF
  244.  
  245.  ENDIF
  246.  
  247.  
  248.  
  249.  case VT of
  250.  
  251.  HARTLEY : BEGIN
  252.               FC = max(Var)/min(Var) 
  253.               if (M EQ 0) THEN DF = [R-1,C] else DF = [max(YC),C] 
  254.               P = -1
  255.              END
  256.  BARTLETT: BEGIN
  257.              if M NE 0   THEN v = Total(YC) ELSE BEGIN
  258.                v = R*C-C
  259.                YC = replicate (R-1,C)
  260.                CN  = replicate(1.0,C)
  261.              ENDELSE
  262.              t1 = N_Elements(Var) -1.0
  263.              CC = 1. + 1.0/(3*t1) *  $
  264.                                        (Total(1./YC(CN)) - 1./v)
  265.              s2 = Total(YC*Var)/v
  266.              FC = 1/CC * (v * Alog(s2) - Total(YC* Alog(Var)))
  267.              DF = t1
  268.              P  = 1. - chi_sqr1(FC,DF)
  269.  
  270.              END
  271.  
  272.  BOX:     BEGIN
  273.             if n_elements(gn) EQ 0 THEN gn = 2
  274.  
  275.             if gn gt R then BEGIN
  276.                print, "equal_variance- group size is too large"
  277.                return
  278.             ENDIF
  279.  
  280.             group = fix(gn * randomu(seed,C,R))    
  281.                                         ; break columns into groups
  282.             if (M NE 0) THEN group(NotGood) = -1     
  283.             Var =  Fltarr (C,gn)
  284.             for i = 0L,gn-1 DO BEGIN    
  285.                                ;compute group variances for columns
  286.               tempgroup = X
  287.               ngri = where ( group ne i)
  288.               tempgroup(ngri) = 0
  289.               gri = where(group eq i)
  290.               tempgroup (gri) =1
  291.               ysize = tempgroup # Replicate(1,R)
  292.               ysize1 = ysize
  293.               y = where(ysize ne 0,count)
  294.               if count GT 0 then BEGIN
  295.                 ysize(y) = 1./ysize(y) 
  296.                 ysize1(y) = ysize1(y) -1
  297.                ENDIF
  298.  
  299.                 tempgroup(gri) = X(gri)
  300.                 Var1 = (tempgroup*tempgroup) # Replicate (1.,R) -$
  301.                          (tempgroup # Replicate(1.,R))^2 * ysize
  302.  
  303.                 y = where(ysize1 ne 0,count)
  304.                 if count NE 0 THEN ysize1(y) = 1/ysize1(y)
  305.                 Var1 = Var1 * ysize1
  306.                 Var(*,i) = Var1        
  307.                                ;append variances for ith group
  308.             ENDFOR          
  309.             
  310.             DF = [C-1,C*gn-C]
  311.            not0 = where (var ne 0,count)
  312.            if count NE 0 THEN var(not0) = alog(var(not0))
  313.            FC = 1
  314.            anova,var,FCTest = FC, ONE_WAY=1  ,/No_Printout
  315.            P = 1 - f_test1(FC, C-1, C*gn-C)
  316.             END
  317.  
  318. LEVENE :  BEGIN
  319.             if ( M EQ 0) THEN BEGIN
  320.                 YC = Replicate (R-1,C)
  321.                 ysize = 1./YC
  322.             ENDIF  else begin
  323.                 y = where(YC NE 0,count)
  324.                 ysize = YC 
  325.                if count GT 0 then ysize(y) = 1./ysize(y)
  326.                 
  327.              ENDELSE
  328.        
  329.             Dev= ABS(X - ((X # Replicate (1.,R))*ysize) #     $
  330.                                            Replicate(1.,R))
  331.             if (M NE 0) THEN BEGIN
  332.                Dev(NotGood) = -1
  333.                tot = total(YC)
  334.             ENDIF ELSE tot = R*C   
  335.  
  336.             DF = [C-1,tot-C]
  337.             FC=1
  338.            if ( M NE 0) THEN      $
  339.                anova,Dev,FCTest = FC,      $
  340.             Unequal_ONE_WAY=1,Missing=-1,  /No_Print $
  341.            ELSE  anova,FCTest = FC,Dev,ONE_WAY=1, $
  342.                       /No_Print          
  343.             P = 1 - f_test1(FC, C-1, tot-1)
  344.             END
  345.             
  346.      
  347.  ELSE :
  348.  ENDCASE
  349.  
  350. RETURN
  351. END 
  352.